#include "dyn_multifunctor_multiMaps_henonMap.h"

HenonMap::HenonMap(const double a, const double b, const int p)
    :MultiMaps(p)
{
    setA(a);
    setB(b);
}

colvec HenonMap::iterate(const colvec X, int p)
{
    // If the period is null, we return the input vector itself.
    if (p==0)
        return X;

    // We use auxiliary variables to update the
    // old variables to new ones.
    double x_temp;
    double y_temp;

    // We declare the return vector and set it
    // equal to the input vector.
    colvec answer(2);
    answer = X;

    // We compute the pth iterate of the map.
    for (int i=0;i<p;i++)
    {
        x_temp = 1+answer(1)-a*pow(answer(0),2.0);
        y_temp = b*answer(0);

        answer(0) = x_temp;
        answer(1) = y_temp;
    }
    return answer;
}

mat HenonMap::jacobian(const colvec X, int p)
{
    // We use an auxiliary vector to hold the
    // pth iterate temporarily.
    colvec temp;

    // We declare the return matrix and an auxiliary
    // matrix to compute the matrix product.
    mat jac(2,2), t_jac(2,2);

    for (int i=0;i<p;i++)
    {
        temp = iterate(X,i);

        // We fill the matrix.
        t_jac << -2.0*a*temp(0) << 1 << endr
              << b              << 0 << endr;

        if (i==0)
            jac = t_jac;
        else
            jac *= t_jac;
    }
    return jac;
}
